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We report self- consistent quasiclassical calculations of 
spontaneous currents and magnetic moments in finite size un- 
conventional superconducting systems, namely: (i) in isolated 
d-wave superconductor islands where, in addition to the dom- 
inant order parameter (with a d x i_ y i symmetry), a subdom- 
inant order parameter of s or d xy symmetry is added; (ii) in 
grain boundary junctions between two arbitrarily oriented d- 
wave superconductors, and between a d-wave and an s-wave 
superconductor. We show that the profile of the spontaneous 
current density and the magnetic field distribution depend on 
the time-reversal symmetry breaking properties of the system. 
For the d x i_ y i -\-id xy state, vortices appear near the edges of 
the finite size systems. We associate these vortices with the 
chiral nature of the mixed order parameter. The method de- 
veloped here is quite general, and can be used for predicting 
properties of any finite size superconducting system. 

The possibility of the existence of mixed order pa- 
rameter symmetries in unconventional superconductors 
has been suggested and recently investigated both the- 
oretically [1-8] and experimentally [9-14]. A signifi- 
cant feature of the mixed symmetry states is that they 
may produce spontaneous currents and magnetic mo- 
ments which can be measured using appropriate exper- 
imental techniques. High T c cuprates present an espe- 
cially interesting class of materials for studying mixed 
symmetry states. In these systems, recently developed 
technology [15] provides an opportunity to fabricate dif- 
ferent structures with controllable characteristics, such 
as 7r-junctions [1,16] (junctions with equilibrium phase 
difference equal to n), submicron size ^-junctions [17] 
(with an equilibrium phase difference </>o which is differ- 
ent from or 7T [5]), tt-SQUIDs [1,18] and tt/2-SQUIDs 
[19], and superconducting qubit prototypes [20]. There- 
fore, a quantitative prediction of properties of such finite 
size unconventional superconducting systems becomes 
extremely important for both experiments and technol- 
ogy. 

The mixed symmetry states can be realized near a sur- 
face of a d-wave superconductor or at a grain boundary 
junction between two different superconductors. In the 
grain boundary junction case, the symmetry of the mixed 
state is dictated by the symmetries of the order parame- 
ters on both sides of the junction (proximity effect) [1-5]. 
On the other hand, at a surface, a subdominant order pa- 
rameter may appear due to an attractive interaction po- 



tential in the corresponding channel. This may involve a 
second order phase transition below the superconducting 
transition temperature [5]. The mixed symmetry states 
that have been suggested for d-wave superconductors are 
d x 2_ y 2 +id xy and d x i_ y i +is [6] (we denote them d + id f 
and d+is, respectively). Which of these states is realized 
near the surface of a d-wave superconductor is still an 
open question. An identification of this state would pro- 
vide useful information about microscopic interactions in 
the superconductor. The d + id' and d + is states are in- 
herently different in the sense that the Cooper pairs in 
the former case have intrinsic magnetic moment, while in 
the latter, they do not [21,5]. Therefore, it is natural to 
expect differences in some measurable quantities, such as 
spontaneous current and magnetic field. 

In this article, we perform self-consistent calculations 
of spontaneous current densities and magnetic field distri- 
butions for finite size systems with different mixed sym- 
metries. We show that these experimentally observable 
characteristics are very dependent on the time-reversal 
symmetry properties of the system. For the d + id' state, 
we find that vortices appear, and we connect their ap- 
pearance with the chiral nature of the order parameter 
in this state. 

Most general approaches to calculation of currents in 
superconducting systems are based on Gorkov equations 
for Green's functions of the superconductor. It is widely 
accepted that this "mean field" approach is valid on a 
phenomenological level [22], independent of further de- 
velopments in the first principles' theory of high T c su- 
perconductivity. In a quasiclassical limit, these equations 
(which are also called Eilenberger equations, see Ref. 
[23]) can be solved by integrating along the classical tra- 
jectories of the quasiparticles (with boundary conditions 
at infinity). In finite size systems, where the trajectories 
undergo multiple reflections from surfaces and interfaces, 
the problem of definition of the boundary conditions be- 
comes complicated. We demonstrate how a modification 
of the Eilenberger formalism based on the Schopohl-Maki 
transformation [24] can be used for stable numerical cal- 
culations in finite size two-dimensional (2D) systems. 

The Eilenberger equations for quasiclassical Green's 
functions can be written as 



v F • V# + [u? 3 + A,g]= 0, g 2 = T, 
where uj is the Matsubara frequency and 
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The matrix Green's function g and the superconducting 
order parameter A are both functions of the Fermi ve- 
locity vp and the position r. A is determined by the 
self-consistency equation, which in the two-dimensional 
(2D) case can be written as 



A(e)=27rN(0)Tj2{Vee>f(e')) 6 



(2) 



oo>0 



where is the angle between vp and the x-axis, Vqq> is 
the interaction potential, N(0) is the density of states at 
the Fermi surface, and {...) e represents averaging over 0. 
We also have A"*" = A*, which holds for any singlet order 
parameter symmetry. Generally, it is possible to obtain 
an order parameter which is a mixture of several terms 
with different symmetries, e.g., A = A x 2_ y 2 + A xy + A s , 
where A x i_ y i = Aicos2#, A xy — A2sin2#, and A s are 
the dominant d x 2_ y 2 component, the subdominant d xyi 
and the s components of the order parameter, respec- 
tively. The corresponding interaction potential, 



Voo' = V d i cos 2(9 cos 20' + V d2 sin 2(9 sin 20' + V s , 



(3) 



should be substituted in the self-consistency equation (2) 
in order to assure a self-consistency of the solution [5]. 
The current density j(r) is given by 



j = -47rie7V(0)T^(v^) e 



(4) 



uj>0 



For numerical calculations, it is convenient to param- 
eterize the quasiclassical Green's functions via [24] 
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Functions a and b satisfy two independent, but nonlinear, 
(Riccati) equations: 



v F • Va = A - A*a 2 - 2ua, 
-v F • Vb = A* - Ab 2 - 2uob. 



(6) 



From these equations it follows that a(— vp) = b*(vp) 
and b(— vp) = a*(vp). The solutions of these equations 
can be interpreted as trajectories of quasiparticles. One 
should integrate these equations along all possible trajec- 
tories and perform the summation over the trajectories 
to calculate the current. To find a and b along the tra- 
jectories, one needs to use boundary conditions at the 
ends of the trajectories. In infinite systems, one usually 
assumes that all trajectories go deep into the bulk of the 
superconductor, i.e., one can use the bulk solutions 



a± 



b± 



(7) 



with Q = ^u 2 + |A| 2 , as the boundary conditions at in- 
finity. In the case of finite size systems where multiple 



reflections from surfaces and interfaces are possible, sta- 
bility of the numerical procedure for calculating the cur- 
rent is not obvious because of complications with choos- 
ing a proper boundary condition for a given trajectory. 
Nevertheless, we show that a stable numerical procedure 
can still be developed. 

Numerical integration for a (b) in Eq. (6) should be 
taken in the direction of vp (-vp) to ensure stability. 
When A is a constant, the solution for a can be written 
analytically, 



af ~ a+ l + 4f(a;-a_)e^sinhftT 

ft ( ' a,i — a + 
— I 



(8) 
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where and aj are the values of a at the initial (r^) 
and final (rj) points of the trajectory, and r = |rj — 
Yi\/vp is the migration time between the initial and final 
points. It is clear that the solution for a relaxes to the 
bulk value a + at the distance L = vp/2Q which is of 
the order of the coherence length £o- I n other words, 
when the quasiparticle moves away from the initial point 
at a distance of a few £o's, any information about the 
initial point en is lost. The same argument is valid for 
the function b. 

Let us now consider a system with a restricted geom- 
etry. After integrating over a few £o's (considering the 
reflections), af becomes almost independent of a^. This 
solution corresponds to a simple exponential relaxation 
of the functions a and b to their local "steady-state" val- 
ues (which can be different from the bulk values) defined 
by the local values of the order parameter. Therefore, to 
find A at a given point, one does not need the values of 
A at distances larger than several £o along the trajectory. 
Such a "relaxational" property of the a and b functions 
significantly simplifies the numerical solution of the self- 
consistent 2D problem. All the "multiple-reflection his- 
tory" of a trajectory becomes a moot point, and for prac- 
tical calculations, this trajectory can be cut at a distance 
of a few £o's from a given point. We used this observa- 
tion in our self-consistent calculations. For the "cutting" 
distance, we choose 10£o-20£o- We set the bulk values 
of a and b (a+ and 6+) calculated at the initial point of 
the "truncated" trajectory as the boundary conditions 
(such a choice does not affect the final results because 
the system has no memory) and integrate along the tra- 
jectory until we get to the point where we calculate the 
current. We found that the results are very stable and 
do not depend on the value of the "cutting" distance. 

To compare the spontaneous currents and magnetic 
fields in the d + is and d + id' states, we performed self- 
consistent calculations of the order parameter in small 
size triangular regions of clean d-wave superconductor 
with different subdominant order parameters. We as- 
sume specular boundaries with the length of the longer 
edge of the triangles being 30£o- Figs. 1 and 2 show 
the spontaneous current and magnetic field distributions 
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FIG. 1. Spontaneous current (left) and magnetic field 
(right) in a finite size triangular-shaped region of a d-wave 
superconductor with a subdominant s-wave, at the temper- 
ature T — 0.2T CS — 0.1T C . The orientation of the dominant 
d x 2_ y 2 order parameter is shown. 

in systems with the two different symmetries. We im- 
posed the desired subdominant symmetry by introduc- 
ing an additional attractive interaction potential in the 
corresponding channel [Vd2 or V s in Eq. (3)]. The sub- 
dominant critical temperatures are taken to be 0.5T C in 
both cases. The orientation of the dominant order pa- 
rameter makes a 45° angle with respect to the longer 
(right) edge of the triangle (see the figures). As a result, 
only this boundary of the triangle is pair breaking; the 
quasiparticles face a different sign of the dominant or- 
der parameter after reflection from the surface [5]. The 
suppression of the main order parameter at the surface 
allows the subdominant order parameter to appear near 
the surface. 

The spontaneous current and magnetic field distribu- 
tions are evidently different in Figs. 1 and 2. In the first 
case, the current makes a counterclockwise loop through- 
out the region, while in the second case, the direction of 
rotation of the current changes throughout the triangle. 
Especially, appearance of the two vortices in Fig. 2 re- 
flects the chiral nature of the d + id! symmetry of the 
order parameter [5]. There is no phase winding around 
these vortices and therefore no flux quantization is ex- 
pected. In fact, the flux trapped in these vortices is 
much smaller than a flux quantum. The magnitude of 
the magnetic field at the maximum positions is of the or- 
der of 10 _4 -10 _3 G in both cases. This value agrees with 
the magnitude of the magnetic field reported in Ref. [12]. 
In the d + id' case, the magnetic field is strongly peaked 
at the vortices. The magnetic field in other spatial points 
is noticeable only in the vicinity of the pair-breaking edge 
(within a few coherence lengths) and is almost one order 
of magnitude weaker than the field at the vortex peaks. 
This is a result of the superscreening effect which hap- 
pens when the spontaneous current is due to the intrinsic 
angular moment of the order parameter [4,5]. The total 
flux generated by this magnetic field is very small and 



FIG. 2. The results of self- consistent calculations similar 
to those of Fig. 1 but with a subdominant d xy -wa,ve order 
parameter, and at the temperature T — 0.2T C 2 = 0.1T C . 

difficult to measure. This is in agreement with the mag- 
nitude of the flux, calculated using a different technique, 
in Ref. [25]. In the d + is case, on the other hand, the 
peak is spread over a wider region. In fact, in larger sys- 
tems, the size of this region could be of the order of the 
penetration depth. This is because in the d + is case, the 
superscreening effect is absent and the only mechanism 
to return the current is the Meissner effect, which hap- 
pens over the length scale of a penetration depth. The 
flux generated by this magnetic field can be large and, 
therefore, measurable. This suggests that the flux mea- 
sured in Ref. [12] can be generated by a d + is symmetry 
breaking state and not by a d + id' one. 

The spontaneous magnetic moments of mesoscopic is- 
lands with d + is or d + id' symmetry can be in principle 
used as qubits. Since the "up" and "down" states are de- 
generate eigenstates of the effective Hamiltonian, tunnel- 
ing between them becomes only possible in the presence 
of external magnetic field, field gradient, etc., depending 
on the multipole structure of the spontaneous current 
distribution. In order to determine the tunneling ampli- 
tude, a three-dimensional picture should be considered, 
which is beyond the scope of this publication. 

We also have calculated the spontaneous current and 
magnetic field distributions in grain boundary junctions. 
The geometry of these systems consists of a square of 
size = 30£o x 30£o> divided into two equal parts sepa- 
rated by a grain boundary junction. The right half of 
the square is a clean d-wave superconductor with a 45° 
orientation of the order parameter with respect to the 
boundaries. The left half, on the other hand, can be ei- 
ther a clean s-wave superconductor (s-d junction, Fig. 3), 
or a clean d-wave superconductor with a 0°-orientation 
(d-d junction, Fig. 4). The grain boundary is taken to be 
perfectly transparent with no roughness or faceting [26]. 
Although this choice of grain boundary does not exactly 
correspond to the reality, study of it may provide use- 
ful information for real systems. To be able to compare 
the two systems, we assumed that the s-wave and d-wave 
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superconductors have the same transition temperatures. 
We did not introduce any subdominant order parameter 
here. However, we introduced an additional phase differ- 
ence of Acj) = 7r/2 between the two sides. Such a choice 
corresponds to the equilibrium phase difference of the 
junction at which the total current passing through the 
junction is zero [5]. Calculations are done at T = 0.1T C . 

The results of spontaneous current density and mag- 
netic field distributions are shown in Fig. 3 for the s-d 
junction, and Fig. 4 for the d-d junction. Note that the 
current distribution is not symmetric with respect to the 
grain boundary. For the s-d junction of Fig. 3, the cur- 
rent has a maximum at the grain boundary and returns 
through the bulk of the superconductors on both sides. 
If the system is large, the size of the region where the 
spontaneous current is non-zero, should be of order of the 
penetration depth. For the d-d junction (Fig. 4), on the 
other hand, the current changes the direction just within 
a few coherence length from the boundary, again due to 
the superscreening effect [4,5]. Notice also that near the 
edges of the system, on the left side of the grain boundary 
junction (0°-orientation of the order parameter), the cur- 
rent returns along the diagonal, whereas on the right side 
(45°-orientation) it forms two small vortices and antivor- 
tices. These vortices, which are absent in the s-d case, 
again reflect the chiral nature of the d + id' symmetry; 
although the subdominant order parameter is absent in 
this case, the correlation functions convey this symmetry 
near the grain boundary junction due to the proximity 
effect [5]. The magnetic field distributions, displayed in 
Figs. 3 and 4, is of the order of 10 _4 G at the positions of 
the maxima, in both s-d and d-d junctions. In the case of 
d-d junction (Fig. 4), magnetic field is peaked at the lo- 
cation of vortices. Away from the vortices, the magnetic 
field is localized near the grain boundary junction and al- 
most one order of magnitude smaller. The flux generated 
by such a magnetic field is very small. Thus, observation 
of large flux as in Ref. [10] can not be associated with 
such an effect. In real systems, due to the faceting effect 
[26], some n- loops may exist along the grain boundary 
junction. They can produce large fluxes, of the same 
order as observed in Ref. [10]. 

It is important to emphasize here that the existence 
of the spontaneous current at the grain boundary junc- 
tions does not depend on the presence of a subdominant 
order parameter (unlike in the previous case, unless if 
we consider other effects [27]). Thus, experiments simi- 
lar to that of Ref. [11] (that only probes the mixing of 
the symmetry of the order parameter) can not exclude 
the possibility of the existence of spontaneous current. 
Indeed, addition of a subdominant order parameter will 
suppress the spontaneous current at the boundary (see 
Refs. [4,5]). 

One should note that in the presence of magnetic field, 
the procedure described here is not valid, because a path 
dependent phase will be accumulated to a and b func- 
tions, and the relaxation mechanism along the trajectory 
may no longer hold. This does not mean that we cannot 
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FIG. 3. Spontaneous current density (left) and magnetic 
field (right) distributions for an s-d junction. The grain 
boundary is a vertical line located in the middle. 




FIG. 4. Spontaneous current density (left) and magnetic 
field (right) distributions for a grain boundary junction be- 
tween two d-wave superconductors. 



use our approach to calculate currents and fields in realis- 
tic systems without calculating the magnetic field due to 
spontaneous currents self-consistently. Indeed, it is not 
difficult to show that the corresponding Doppler shift is 
rather small. We saw that for systems with size L ~ 
30£o - 5 x 10 4 A, spontaneous magnetic field is always 
smaller than 10 _3 G. Therefore, the induced superfluid 
momentum is at most p s ~ (e/c)HL ~ 10 _26 g • cm/s. If 
we take Tivp — 1 eV-A [28] for the Fermi velocity in high- 
T c cuprates, we obtain p s vp ~ 0.1 fieV ~ 1 mK, which is 
negligible compared to other energy scales in the problem 
[29]. However, these restrictions may be important if the 
system is placed in a strong external magnetic field. 

To summarize, we described self-consistent calcula- 
tions of some equilibrium properties in finite size su- 
perconducting systems. We calculated distributions of 
spontaneous current and magnetic field in different small 
samples in the presence of mixed order parameter sym- 
metries. The nature of the mixed symmetry state in all 
cases affects the shape of the current and magnetic field 
distributions. In particular, the chiral nature of the d+id' 
states exhibits itself through the appearance of vortices 
close to the edges of the system. The vortices are ab- 
sent in the d + is cases. The method described here is 
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quite general and can be applied to any 2D geometry with 
proper boundary conditions as long as external magnetic 
field is not present. The self-generated magnetic field by 
the spontaneous currents, however, is usually very small, 
i.e., it can be neglected in most of the practically impor- 
tant cases. 

Many discussions with A. Maassen van den Brink, 
G. Rose, and A. Yu. Smirnov are greatly acknowledged. 
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